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ABSTRACT 


The  generalized  inverse  is  of  increasing  importance  for  estimation 
and  optimization  in  modern  systems  theory,  because  it  both  simplifies 
many  problems  and  reveals  underlying  structures  of  theoretical  impor- 
tance. Optimization,  using  the  gradient  of  the  trace  of  products  of 
matrix  valued  functions,  including  the  generalized  inverse,  are  pre- 
sented in  a novel  state-space  setting.  A number  of  functions  <jf  such 
products  of  matrices,  including  general  formulas,  are  derived  for  the 
first  time  to  our  knowledge.  Tables  of  a number  of  them  are  given  with 
some  applications  in  the  fields  of  estimation  and  optimization. 
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NEW  GRADIENT  TECHNIQUES  FOR  TRACES  OF  FUNCTIONS  OF  RECTANGULAR  MATRICES 

AND  THEIR  PSEUDO- INVERSES 


I.  INTRODUCTION.  The  state  space  approach  to  modern  estimation  and  con- 
trol theory  requires  optimal  feedback  or  weighting  matrices  which  mini- 
mize the  trace  of  quadratic  matrix  functions.  In  the  estimation  problem 
the  quadratic  matrix  function  is  the  variance  of  the  state  estimation. 

For  example,  in  the  Kalman  theory  the  best  estimate  requires  that  the 
error  matrix,  P,  be  minimized  with  respect  to  a weighting  matrix  W,  where 
P propagates  via  the  matrix  P.iccati  equation. 

P = AP  + PA^  - WHP  - PH^w”^  + WRw"^ 

The  novel  vector  and  matrix  partitioning  techniques  developed  in  this 
paper,  eliminate  much -of  the  tedium  of  the  analysis  based  on  scalar  sums 
and  products.  In  addition,  they  reveal  a depth  of  structure  not  other- 
wise apparent . 

Section  II  is  devoted  to  deriving  gradient  relations.  Initially, 
considerable  detail  is  exhibited  to  show  necessary  relationships;  after- 
wards, this  detail  is  bypassed.  A table  of  gradients  for  a number  of 
common  important  forms,  principally  linear  and  quadratic  fonris,  is  in- 
cluded at  the  end.  Athans  [1]  presents  some  linear  gradient  forms  for 
the  trace  as  well  as  some  for  determinants,  logs  and  exponentials;  Tracy 
and  Dwyer  [2]  present  additional  gradient  forms  some  of  which  involve  the 
Kronecker  matrix  product;  and  Neudecker  [3]  presents  gradients  involving 
traces,  latent  roots  and  determinants  which  are  related  to  frequently 
occurring  statistical  econometric  problems. 

In  Section  III  a general  formula  is  developed  for  functions  express- 
able  as  a power  series.  Particular  reference  is  made  to  trigonometric 
and  hyperbolic  functions. 

In  Section  IV  applications  for  estimation  and  filtering  are  demonstr- 
ated by  developing  criteria  for  the  continuous  Kalman  filter. 

II.  GRADIENT  DERIVATIONS.  Several  useful  relations  involving  the  grad- 
ient of  some  scalar  valued  functions  with  respect  to  vectors  and  matrices 
are  developed  in  this  section.  The  cases  ar«  restricted  to  linear,  quad- 
ratic and  "cubic"  matrix  terras;  the  latter  implying  the  generalized  in- 
verse. Throughout  the  derivations,  use  will  be  made  of  well  known  results 
of  traces  of  matrices,  namely 

tr(A+B)  = trA  + trB 


I 


(1) 


tr  M = tr  M 

and  the  cyclic  property  of  products.  For  example 

trRST=trTRS=trSTR  (2) 

pxm  mxs,  j,xp  ixp  pxm  mxi  mxS.  ixp  pxm 

A.  LINEAR  CASES.  Consider  the  ix^  matrix  function  of  the  rectangular 
matrix  X of  size  pxm,  that  is 

“"I  ^ (3) 


with  trace 


tr  tr  AXB  (jjote  this  Z is  not  the  matrix  size.) 


the  mxm  matrix  Lr 


BAX 

mx)l  Jixp  pxm 


with  trace 


tr  L^  = tr 


and  the  pxp  matrix 


Lg  = X B A 

nxn  P*^"'  ^^P 


and  by  Equation  (2) 


!,  = tr  AXB  = tr  BAX  = tr  XBA 


The  differential  of  the  trace  is 


dfc  = d tr  (AXB)  » d tr  (BAX)  = d tr  (XBA) 


de  = tr  A dX  B = tr  BAdX 
Hxp  pxm  mxi  mxm 

= tr  (dXBA) 
pxp 


L 


The  differential  of  the  ra>«n  matrix  can  be  written  as 


dl2  = BA  dX  - — Y dX 

mxp  pxra  3X  pxin 
raxm  ^ 

mxp 

and  the  differential  of  the  pxp  matrix  can  be  written  as 

dL  = dX  BA  = dX 

pxm  mxp  pxm  x" 
mxp 

The  differential  of  the  traces  are 


di  = tr 

■"3)1 

T 

dX 

= tr 

dX 

3)1 

737’ 

3X 

mxp 

pxm 

3X 

_pxm 

, 

mxp 

and  by  Equation  (lO)  and  Equation  (13) 
3)1 

= B A 
3X  mxJZ,  Jlxp 

or 


tr  (AXB) 


tr  (BAX)  = (XBA)=  B A 
3X  mxm  3X  pxp  nixil  g,xp 


3X  ILxl 

For  the  special  case  when  the  Li  are  scalers,  i.e.;  A and 
^ tr  ^)  a X b (^  = b(i^^)a 


3X 

3tr 

3X 


pxm 


j^b(^^)a  ^ J" 


and 


(11) 

(12) 

(13) 

(14) 

(15) 

are  vectors 

(16) 

(17) 


3 


[pL  = b(^p)a 


By  Equation  (3)  for  m=Jl  and  B=Im 

^1  A X \ = A 
3X  \ mxp  pxm ) mxp 
mxp 

By  Equation  (3)  for  t=p  and  A=I 


X B ) = 
3X  \ pxm  mxp  / 


In  the  above  equation  if  ra=l 


= 4p)b 


Since  the  trace  of  the  dyadic  product  is  the  inner  product 
tr  (3^^)  = = I 


one  has 


Kk  = <1^  = 4 


where  the  matrix 


2-1  P J 


when  the  coordinates  are  all  independent. 
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One  can  verify  all  of  the  previous  results  via  the  tedious  process 
of  pctrtitioning . For  example  consider  Equation  (20) 

L = X B (25) 

pxp  pxm  mxp 

and  partition  X into  its  column  vectors  and  B into  its  row  vectors  to 
obtain 


XB 


= [x(p>^,...x(p>J 


^^)b 


i=l  ^ 


U“<p)bJ 


the  trace  of  Equation  (26)  is 

tr  XB  = tr  I ( ^2^^  ) 


+ • • • + tr 


(«») 


or 


a = tr  XB  = ^4x>^  + + • ♦ • + /bX 


= + ^2  ^ • • • + 


The  differential  of  t is 


dH  = dt,  + dt„  + • • • + dt 
12  m 


(26) 


(27) 

(28) 

(29) 

(30) 


where  each  of  the  is  a function  of  the  vector  or  by  Equation  (l3) 


.31. 


T— d3^. 

^ 3x  / 1 

Repackaging  Equation  (31) 


(31) 


/3t.  ydl 

‘ \3ir 

1 m 


(32) 
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n 


By  Equation  (31)  in  Equation  (28) 


(33) 


Consider  the  product  of  the  matrix  of  gradient  vectors  with  the 
matrix  of  dX  vectors  that  is 


mxp 


Clearly 


(34) 


(35) 


(36) 


(37) 


The  other  results  can  be  shown  by  similar  partitioning  and  can  be 
found  in  detail  in  reference  [4]. 
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S’  QUADRATIC  CASES.  Consider  the  matrix  product 

..T 


C X 


B 


or 


ixi 


= Tax  irc  x’’  b1 

jix^  J 


(38) 


(39) 


Form  the  new  matrix  of  size  mxm  from  Equation  (39)  cyclically  permuted 
via  parentheses 


Q2  = []c  x^  bJ  ax 

mxt  txm 


”‘2 
mxiti 


(40) 


Form  a th^d  matrix  of  size  pxp  by  cyclically  permuting  the  A matrix  of 
Scjus'tion  (39) 

Q3  = X C x'^  B A (41) 

pxp  P*"*  '"^P  P’^^  ^'‘P 

The  differential  of  Equation  (40)  is 

dQ2  = C dX^BAX  + Cx'^BAdX  (42) 

mxm 

The  first  matrix  on  the  right  of  Equation  (42)  under  the  trace 
transposition  and  permutation  rules  becomes 

tr  (Cdx'^BAX)  = tr  C^A^B^dX  (43) 

Using  Equation  (43)  in  the  trace  of  Equation  (42)  yields 

tr  dQ2  = tr  (C^x\^B^+CX^BA)  dX  (44) 

mxm 


The  gradient  factors  of  dQ^  can  be  taken  to  be 
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If 


BA  = I 

then  Equation  (51)  becomes 

= 2 x^ 

3X  3X  mxp 

For  the  special  case  when  Equation  (49)  is  a scalar,  one  obtains 

[<aXCx'^b>]  = C V^b  + CX'^^4 
3 X 

and  for  C=I 


(52) 


(53) 


(54) 


^[<aXX^b>=  [|><b»b><.3 


3X  mxp 

For  p=l  by  Equation  (51) 


(55) 


3 

— 1 
3X 


<^BA3^  = j 


(56) 


For  C=I  in  Equation  (54) 


3tr 

~T 

3X 


(57) 


As  a final  quadratic  case  consider  the  (full  rank)  grammian-matrix 


T 

X X = G 


and  Its  inverse 


(x'^'x)"^  = g'^ 


(58) 


(59) 


where 


G "G  = I 


(60) 
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The  differential  of  Equation  (60)  is 

dG'^G  + G"^dG  = 0 (6i) 

or 

dG"^  = -G"^dGG”^  (62) 

= -G'^CdX^X+x'^dX]  g"^  (63) 

or 

-1  -1  T -1  -IT  -1 

dG  =-G  dX  XG  - G X dXG  (64) 

The  generalized  inverse,  for  the  full-rank  case,  is  defined  as 

X*  = (X^X)"^  x"^  = g‘V  (05) 

mxp 

and  its  transpose  is 

X*"^  = X G-^  (66) 

pxm  pxm  mxm 

Using  Equation  (65)  and  Equation  (66)  in  Equation  (64) 

dG"^  = - G~^  dX^X*"^  - X*dXG"^  (67) 

The  trace  of  Equation  (67)  is 

tr  dG"^  = - tr  (dx’^X*’^G'^>  - tr  (G~^X*dX) 

Using  the  trace  transpose  property  on  the  first  right  hand  side  term 
of  Equation  (68) 

tr  dG“^  = - tr  (G"^X*dX+G'^X*dX)  (69) 

= - 2 tr(G"^X*dX) 
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or 


(x'^X)'^  = - 2 (X^X)"^X* 
8X  mxm 


Using  Equation  (65)  in  Equation  (71) 


^ (x'^x)-i  = 


2 (X^X)"^  X^ 
mxm  mxp 


C.  "CUBIC"  CASE.  The  following  cases  do  not  involve  cubics  but  the 
matrix  X and  X appear  three  times  in  the  generalized  inverse  relation 
(full  rank) 

X*  = (X^X)‘^  x'^  (73) 

mxp  mxm  mxp 

Consider  the  linear  form  in  X*  that  is 

Q = X*  B (74) 

mxm  mxp  pxm 


dQ  = dX*B 


dQ  = d C(x'^X)“V]  B 


= d(X^X)~^  X^B  + (X^X)"^  dX^B 


by  Equation  (63)  in  Equation  (76) 


dQ  = - G'^  [dx'^X+x'^dX]  G“^X^B  + G'^dx’^^B 


or  by  Equation  (73)  in  Equation  (77) 


dQ  = - G"^[dX^X+x’'dX]X*B  -^G'^dX^B 
Using  the  relation  for  the  projector 

T 9 

X X*  = P = P = P 

pxm  mxp  pxp 

in  Equation  (78) 

dQ  = - G"^  dX^  PB  - X*dXX*B  + G“^  dX^B 
or 

dQ  = g’^  dX^  (I-P)  B - X*  dXX*B 

The  orthogonal  compliment  projector  is 

-T  -2 

P = I - P = P = P 

Using  Equation  (82)  and  Equation  (81) 

dQ  = G'^  dx’^PB  - X*dXX*B 

The  trace  relations  yield 

tr  dQ  = + tr(dX^PBG”^)  - tr  (X*BX*dX) 
and  transpose-wise 

tr  dQ  = tr  ( G‘V?dX)  - tr  X*BX*dX 
or 

tr  dQ  = tr  [(  G'^B^P-X*BX*)dX] 

Expressing  dQ  as 

dQ  =ii£^dx 
mxm  3X  pxm 
mxp 


(78) 

(79) 

(80) 
(81) 

(82) 

(83) 

(84) 

(85) 

(86) 

(87) 
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one  obtains 


(X*B)  = 

G‘^  B^P  - X*  BX* 

9X  ram 

(88) 


D.  POWERS  OF  (X*B)  . The  following  cases  are  of  higher  degree  than 
third.  Consider  n=2  then  the  Q of  Equation  (74)  has  square 

= X*BX*B  (89) 

and 


dQ^  = dX*BX*B  t X*BdX*B  (90) 

By  Equation  (75) 

dQ  = dX*B  (91) 


hence  Equation  (91)  in  Equation  (90) 

2 

dQ  = dQX*B  + X*BdQ  (92) 

By  Equation  (89)  let 

q2  = tr  (93) 

then 

2 ^^2 

dQ  = — - dX  = dQQ  + QdQ  (94) 

3X^ 

The  trace  of  Equation  (92)  is 

tr  dQ^  = tr  2X*BdQ  (95) 

By  Equation  (94)  and  Equation  (95) 
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r 


1 


2 

^ tr  dQ  = tr  Q dQ  = tr  X*BdQ 
Using  dQ  given  by  Equation  (83)  in  Equation  (96) 


tr  X*BdQ  = tr  {X*B  (G"^dX^PB-X*dXX*B)} 

= tr  X*BG  ^dX^pB  - tr  X*BX*dXX*B 
= tr  dX^  (pBX’i'BG"^)  - tr  X*BX*BX*dX 
= tr  (G  ^B^X*'^B'^f  )dX-tr  Q^X*dX 
Js  tr  dQ^  = tr  ((g'^b\*Vp-Q^X*)  dX} 


or 

9q 

*2  —1  T O 

h — s-  = G Q^B  P - Q X* 

ax 

That  is, 

^^2  -ITT  9 

— = 2 (G  Q B P-Q  X*) 

ax 

Let  n=3  and 

= QQQ 

with 

tr  » q3 

The  differential  of  Equation  (103)  is 

3 2 2 

dQ  = dQQ  + QdQQ  + Q dQ 


(96) 

(97) 

(98) 

(99) 

(100) 

(101) 

(102) 

(103) 

(104) 

(105) 

(106) 
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and  the  trace  of  Equation  (105)  is 


3 

tr  dQ 

= 3 tr  Q^dQ 

(107) 

Equation  (83)  in  Equation  (107) 

3 

tr  dQ 

= 3 tr  [G~^dx'^PB-X*dXQ] 

(108) 

Using  the  commuting  and  permuting  properties 

of  trace  one  obtains 

^=3 

V ' 

[g"^  (Q^)^  b^?  - 

(109) 

mxp 

For  n=l,2,  and  3 by  Equation  (88),  Equation 

(101)  and  Equation  (109), 

-IT' 

—f  = G 
3X 

: BP  - QX* 

(110) 

~-l^T„T'  2 

, 

T ^ 

3X 

G Q B P - Q X*  1 

(111) 

T,2  T'  3 ri 

3X^ 

G (Q  ) BP  - Q X* 

(112) 

This  gives  the  inductive  step,  therefore  for 

Q = X* 

B 

(113) 

mxm  mxp  pxm 

-1  T -1 

G ^ = (x'x)  ^ 

(114) 

mxm  mxm 

X*  = g” 

mxp 

(115) 
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P = XX* 

pxp 

p = I - p 

Si  “ tr  q";  n = 


we  have 


3X 


~ = n [G‘^(q’’)" 


" 1 T 

B P-q"x*3 


(116) 

(117) 

(118) 


(119) 


A number  of  special  cases  follow 
square  and  full  rank,  then 


from  Equation  (119). 


Suppose  X is 


and 


Then 


and 


X*  = X"^ 
mxm 


Q = X'^B 


P = XX‘^  = I 


P = 0 


hence 


^ (X-S)'^  = - n (X-^B>"  X-1 


3X 


Suppose  now  that  B=I  then 

(X-^)"  = - n (X'S"  X"^ 


3X 


(120) 


(121) 


(122) 


(123) 


(124) 


(125) 
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^ (X-^"  = 


E.  POWERS  OF  (XB)  » The  powers  of 
Q"'  = (XB)" 


can  be  obtained  in  a similar  manner.  By  Equation  (20) 
3qi 


— T ' ® 
3X  mxp 


2 2 
Q = (XB)^ 


qo  = tr  Q 


dQ"^  = dQQ  + QdQ 


tr  dQ  =2  tr  (QdQ) 


= 2 tr  QdXB  = 2 tr  (BQdX) 


= 2 BQ 


Repeating  the  arguments  as  before  one  obtains 


Ctr  (XB)"]  = n B (X  B )"'^ 


mxp  pxp 


(135) 


For  X square,  full  rank,  and 


B = I 

Equation  (135)  becomes 

3trX^  ,.n-l 
If  = nX 

3X^ 


F.  TABLES  OF  GRADIEOTS 
!•  Linear  Forms: 

3tr  / . ' 


3tr 

...T 


3tr 

3X^ 

3tr 

T 

3X 

3tr 


3tr 

T 

3X 

3tr 


\ 

pxm 

mxi  / 

i ® 

A 

X ) 

\ mx£ 

Axp 

pxm  / 

/ X 

B 

A ) 

X 

B 

mx4 

Jlxp  / . 

[b(r^^)a 

X ] 

pxm 

[ X 
Pxm 

b(i^^p)a]  1 

[ A 

X ] 

= A 

mxp 

pxm 

mxp 

C X 

B ] 

= B 

X 

3 

mxp 

mxp 

} = BA 


^(x(^4)b)  = <p)b 
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1 


(136) 


^ 


2.  Quadratic  Forms: 


3tr 


3tr 

„,.T 


3tr 

...T 


3tr 


3X 

3tr 


3X 


[ A 

X 

C 

T 

X 

B ] 

txp 

pxm 

mxm 

raxp 

pxA 

[ c 

T 

X 

b 

A 

X J 

mxm 

mxp 

pxi 

£xp 

pxm 

[ X 

C 

X^ 

B 

A ] 

pxm 

mxm 

mxp 

pxn 

Zxp 

[Xx'^^BA] 

CAXX"^ 

B] 

' = 

)^(A^ 

b'^+ba) 

[x’^bax] 

T T T T T 
= C X A b'  + Cx'bA 


3X 


■—  [x'^X]  = 2x’'’ 

3X 


^ «aXx''^^) 

(x'b^x) 

9X 


ax 

3tr  /vT 


= X^  (^<^+:^<a) 


(x'^’x)'^  = - 2(x'^X)“^  x'*' 

3X 
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3. 


"Cubic"  Forms  and  Others:  The  generalized  inverse 


X*  = (X^X)"^  X^  = G'^ 

^ (X*  B ) = (x\)-^  b'^P  - X*BX* 

3X  mxp  pxin 

where  the  projectors  are 
P = I - P 
P = XX* 

G = X^X 

^ (X*B)"  = n [G"^(q’^)““^  b'^P  - Q^*] 

3X 

~ (X“^B)"  = - n (X“S)"  X“^ 

3X* 

- n 

nB(XB)""^ 

^ (X")  = nx"-l 
3X 

III-  S^QME  FUNCTIONS  OF  GENERALIZED  MATRIX  PRODUCTS  INVOLVING  (x  x'^  X*  X*'^^ 
OF  FINfTE  ORDER.  In  this  section  we  wiU  be  concerned  with  matrix  fincUons 
of  the  form 


(X'l)"  = 


3X 


3tr  (XB)  = 

T 

3X^ 


Y - 


•B  -Z  -B  , 
n n n+1 


n>o,  n an  integer  and  the  Bj  scalar  matrices.  Functions  of  Y,  the  argu- 
ment, will  be  examined,  such  as  f(Y)=Y,  and  certain  transcendentals 
identified  below.  The  algorithms  to  be  developed  are  divided  into  three 
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classes,  viz:  those  in  which  the  Zie{X,X'^},  where  the  Zie(X*,X*T)  and 

those  arguments  where  the  can  be  any  of  these  matrices.  It  will  be 
seen  that,  for  each  of  the  functions  (mappings)  listed  (in  each  class), 
the  partial  derivative  of  the  trace  with  respect  to  XT  has  a common  form, 
and  that  the  different  solutions  merely  require  the  substitution  of  a 
specified  matrix  into  this  form.  The  generalized  equations  will  be  speci- 
fied, but  the  formal  proof  will  not  be  given.  These  equations  can  easily 
be  proved  by  induction,  but  only  the  induction  step  is  shown.  It  will  be 
clear  from  the  foms  of  the  first  two  classes  how  the  results  can  be 
combined  to  solve  for  a generalized  mixed  product.  This  will  be  indica- 
ted and  the  general  solution  presented. 


Finally,  at  the  end  is  a table  of  the  partials  of  the  trace  for  some 
of  the  more  commonly  expected  products. 

The  following  set  definitions  will  be  helpful.  Let 

S,  = {B  : the  B are  finite  rectangular  scalar  matrices.} 

A 

S = {X  : the  X are  full  rank  rectangular  variable  matrices. 

XeS^-f>XTeS  } 

X x 

S = {X*  : the  X*  are  the  pseudoinverses  of  the  elements  S .} 

X X 

Sy  = (Y  : the  Y are  square  matrices  of  the  form 
n 


Y = ( n B.X.)  B 

. , 11  n+1 

1=1 

where  n>o  is  an  integer,  B^eSy^,  XjeSx) 


Sy  = {Y  : the  Y are  square  matrices  of  the  form 


Y = ( n B.X*)  B , 
, T 1 1 n+1 

1=1 


where  n^o  is  an  integer, 


S2  = (Z  ; the  Z are  square  matrices  of  the  form 
n 

Z = ( n B.Z.)  B , 

. , 11  n+1 

1=1 

where  n>o  is  an  integer,  Z^  e S VJ  S*> 
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Class  I differentials  will  be  developed  first  in  considerable  detail 
using  the  argument: 

Y = A X B x'^C  X D;  A,  B,  C,  DeS.;  X,  x'^  S ; YeS„  (1) 

A X Y 

followed  by  a more  abbreviated  development  for  Class  II  differentials 
using  the  argument: 

Y = A X*  B X*  C;  A,  B,  CeS  , X*eS  , YeS*  (2) 

A X Y 

Specific  functions  of  Class  III  arguments  will  not  be  developed. 

The  following  auxilliary  definitions  will  be  helpful 

D^4  AdZ  (Bz'^CZD) 

D^4  (AZB)dz'^  (CZD)  (3) 

D34  (AZBZ'^C)  dZD 


where  ZeS^  or  ZeS*.  The  restriction  or  extension  of  the  subscripts  is 
obvious. 

Let  SeS/^.  Then  (since  tr  (A)=tr(AT)) 
tr  (SDj^)  = tr  (BZ^CZDSAdZ) 
tr  (SD^)  = tr  ((AZB)"^  s'^(CZD)'^  dZ) 
tr  (SDg)  = tr  (DSAZBZ^CdZ) 

Let  Fm  4 {±SY,  6'^,  f sinY,  ± cosY , ± sinhY,  ± coshY  : SsSa» 

YcSy  U S$}  and  M be  a mapping  such  that  M : Y-+-M(Y)eFm.  Call  this  function 
in  FM>  Q>  and  in  all  cases,  let  q be  defined  as 
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q 4 tr  (Q) 


Class  I 


Let 


Y = AXBX^CXD 


Q e Fm  = SY, 


SeS 


A 


Then 


dX4dQ  = SCD^+D^+D^) 


(5) 

(6) 


(7) 


where  the  Di  are  given  by  Equation  (3).  Using 
stituting  Equation  (4)  into  Equation  (7) 


trace  properties. 


and  sub- 


tr  (DQ)  = tr  (SD  ) + tr  (SD^)  + tr  (SD  ) 
■‘■2  3 

= tr  [ ( BX'^CXDSA+  ( AXB ) ^s'^(CXD  )^-tDSAXBx’^C  )dX] 

Therefore 


(8) 


■^^  = BX^CXDSA  + (CXD  SAXB)^  + DSAXB  x"^C 
3X 


and  if  S=I 


(9) 


3q  _ 9(trY) 

3?"  ax”^ 


T y 

BX  CXDA  -h  (CXDAXB)  ■t  DAXBx’^C 


(10) 


Let  QeFM=eY.  Then  Q can  be  expanded  in  a series  [5.  6 7 in 
particular,  theorem  4,  p 46,  ref  7]  as  » » » 


I 


Y 


■t 


(11) 


Then 
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dQ  = t + Dg  + ^7  (D^Y+D^Y+D^Y+YD^+YD^+YD  ) 


|r  (D^Y^+D2y2+D3y2+YD^Y+YD2Y+YD3Y+y2d^+y2d  +y^D 


and 


tr(dQ)=tr  [(I+Y+-^Y^+- • • )D  J + tr  [(l+Y+i7Y^+- • • )D  ] 
• J*  2*  2 


+ tr  [(I+Y^i^Y^+•..)D3] 


since  tr(DiY+YDi)=2tr(YDi) , etc.  Therefore 
tr(DQ)=  tr  [e^(D^+D2+D3)] 


2 ■ “3 

(12) 


(13) 


(14) 


Since  e^eSA,  identify  S with  eY  and  substitute  into  Equation  (9).  This 
gives 


= BX^CXDe^A  + (CXDe^AXB)^  + De^AXBx'^C 
3X 


(15) 


We  might  note  that,  due  to  the  permutative  aspects  of  the  trace 

ax 

can  be  written  in  a variety  of  ways.  The  one  chosen  here  appears  most 
"natural"  because  the  argument,  Y,  remains  unchanged.  However,  suppose 
YAAXB  and  Q^e*,  then  the  following  expressions  are  equivalent 


-■X  = Be*’'”  A(.Be>A)  . = BA 

ax 


XBA 


A similar  set  of  equalities  exist  for  each  function  in  Fjj.  However, 
the  representation  is  unique  if  the  argument  remains  fixed,  in  the 
remainder  of  this  section,  only  the  form  which  preserves  SeF(4  will  be 
considered. 
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Now,  had  Q been  eiY=cosY+isinY , the  terms  in  the  expansion  of  dQ  could 
be  collected  into  real  and  imaginary  parts.  It  can  be  easily  seen  that 
these  sums  can  be  written  as  £3] 

tr(dQ)  = - tr  [(sin  YKD^+D^+D^ )]  + i tr  [(cos  YHD^+D^+D^)]  (16) 

Let  qi4tr(cosY)and  q24^(sinY).  Then  if  S is  chosen  as  S=-sinY , Equation 
(9)  becomes 

— Y = - B X CXD  (sin  Y)  A 
3X 

(17) 

- (CXD(sin  Y)  AXB)^  - D(sin  Y)  AXBX^C 
If  S is  chosen  as:  S=cosY,  Equation  (9)  becomes 

^^2  T 

— Y = B X CXD  (cos  Y)  A 
9X 

(18) 

+ (CXD(cos  Y)  AXB)^  + D(cos  Y)  AXBX^C 

Similarly,  if  the  expansions  of  tr(eY)  and  tr(deY)  are  grouped  accord- 
ing to  even  and  odd  powers,  3/3XT  tr(coshY)  and  3/3XT  tr(sinhY)  can  be 
expressed  by  substituting  S=sinhY  and  S=coshY,  respectively,  in  Equation  (9). 

There  is  a more  enlightening  way  the  above  derivatiom can  be  performed. 

Let 

Xj^AAX 

X^(Xb'^)'^  (19) 

X3ACX 

then 


Y = 
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Also,  for  any  i=l,2,3,Y  can  be  expressed  as 


Y = L.X  M. 
1 i 1 


(except  for  a transpose  on  X2)  and 


4 - ' 


"1  = 


4 = "1 


M2  = X3D 


(using  X_) 


(20) 


^ ■ "1^2 


M3  = D 


Then  Equation  (4)  can  be  written  as 

tr(SD^)  = tr  (Mj^SL^dXj^)  = (M^SL^AdX) 

tr(SD2>  = tr  (M2SL2dX2)  = tr  ( (M2SL2B)^dX) 
tr(SD3)  = tr  (M3SL3dX3)  = tr(M3SL3CdX) 


(21) 


Such  expressions  can  be  generalized  as  follows.  Let  a^"!"  or  "T" 
used  as  a superscript.  That  is,  for  any  matrix.  A, 


A = 


A if  a = 1 


(22) 


LA  if  a = T 


Define 
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(23) 


where 


end 


X.  A B.X 
1 “ i 


x{  A (XB^)”^  = By,  B.eS^,  i = l,2,---,n+l 


(24) 


Also,  for  any  i,  i=l,2,---,n,  Y can  be  defined  as 


Y = L.x“m. 

Ill 


where 


Li 


i-1 

n x“ , M . 
j=i  ' ^ 


where  L^AI  and  MnABn+1*  Then  for  SeSa  and  Q=M(Y)eFn, 


tr(dQ)  = tr 


i=l 


(M^SL^B^)”dX 


X 


which  gives 


I 

ax  i=i 


(M.SL.B.)“, 
1 11 


(25) 


(26) 


(27) 


(28) 


Note  that  if  M(Y)e{  TY,e  ^:TeSa}  ^^en  S=M(Y);  if  M(Y)=cosY,  then 
S=-sinY;  if  M(Y)=sinY,  then  S=cosY;  if  M(Y)-sinhY,  then  S=coshY;  and  if 
M(Y)=coshY,  S=sinhY. 


Class  II • The  development  of  a generalized  form  using  the  pseudo- 
inverse parallels  the  previous  derivations  very  closely;  only  the  form  of 
the  final  result  is  changed.  Therefore,  it  will  be  sufficient  for  illus- 
trative purposes  to  use  Y as  defined  in  Equation  (2) 
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Y = AX*BX*C;  A,B,CeS.,  X*eS*,  YeS* 

A X I 

By  assumption  X is  full  rank,  so  Equations  (114)  through  (117)  of  the  pre- 
vious section  apply.  Some  of  these  are 

X*  = (x’^x)'V 
G 4 (x'^x)  -i>X*  = G"V 

X*X  = I 

XX*4P,  a projector 
and  G=g'^,  P-p'^.  Let 

D^4AdX>^BX*C  (3Q) 


D^X*BdX*C 

Tlien  for  SeSa 

tr(SD^)  = tr(BX*CSAdX*) 

^(SD^)  = tr(CSAX*BdX*) 

We  require  dX*.  From  Equation  (29) 

. -IT  -IT 
dX*  = dG  X + G dX  , 

From  equation  (62)  of  the  previous  section 


dG“^  = - G"^dGG“^ 


and  from 
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T 

G = XX 

dG  = dx'^X  + x’^dX 

Thus 

dX*  - - G"^dx'^XG"^x'^  - G’^x'^dXG'^x'^  + G'^dx'^  (32) 

If  TeS^.  then 

tr(TdX*)  = - tr(TG'^dX^XG"^X^)  - tr(TG“^x'^dXG" V) 

+ tr(TG"^dX^) 

= - tr(TG"^dX^XX*)  - tr(TX*dXX*)  + tr  (TG”^dX^) 

(33) 

= - tr(PdXG'^^)  - tr(X*TX*dX)  + tr  (dXG"^^) 

= tr(G"^'^(l-P)dX)  - tr(X*TX*dX) 

Thus,  if  Q=SY,  ScSa,  define 
A BX*CSA 

(34) 

A2  A CSAX*B 

Then 

tr(dQ)  = tr(SD^)  + tr(SD2)  = tr(A^dX*)  + tr(A2dX*) 

(35) 

= tr[G'^(A^+A2)(I-P)dX]  - tr[X*(A^+A2)X*dX: 
or 

= G‘^(Ay+Aj)(I-P)  - X*(A  +A„)X* 

3X  ^ ^ ^ 

We  note  that  if  QeFw,  then  SeF|^  according  to  the  scheme  outlined  for 
Equation  (28). 
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The  above  procedure  can  be  generalized  in  the  same  manner  as  before, 
where  "a"  as  a superscript  has  the  same  meaning  as  in  Equation  (22).  Let 


'^(i"l  rntl 


(36) 


X.  = B.X* 
1 1 


(X*)'^  4 B (X*)'^  = (X*bT)'^,  B.eS,,  and  X*eS* 

1 1 A X 


(37) 


Then,  for  any  i,  Y can  be  defined  as 


Y 


A L.(X. ) M. 

~ 1 1 1 


(38) 


where  Lj^  and  ^ave  the  same  definitions  as  in  Equation  (26)  with  X;. 
It  follows  that 


tr(dQ) 


tr 


“ n ” 

y (M.SL.B.)°*dX* 

Li=l  " " " J 


Now  define 


(39) 


n 

r 4 I (M  SL  B )“  (40) 

i=l 

then 

tr(dQ)  = tr(rdX»'')  = tr[(G"^r^(I-P)-X*rX*)dX]  (41) 

Thus 


= G“^r^(i-p)  - x*rx* 

BX 


(42) 


and  the  same  remarks  apply  as  for  Equation  (28). 
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F 


Class  III.  From  the  foregoing  derivations  it  is  clear  that  if 


n 

n 

i=l 


Y=(  n Z YeS^ 


where  can  be  either  BiX“  or  Bi(X*)“,  then  for  any  i 


Y = L.Z.M. 
Ill 


and,  if  Z ^is  X or  X*,  then 


tr 


(dQ)  = tr  I (M.SL.B. )“dZ' 

Li=i  ^ J 


" "^1 

" "2  -| 

' 1 (M,  SL,  B,  )“dX 

Lk=i  X X 

+ tr 

y (M  SL  B )®dX-’ 
m m m I 

_m=l  J 

(43) 


where  k sums  over  the  ni  terms  involving  dX  and  m sums  over  the  n2  terms 
involving  dX*;  the  Mi,  Li  have  the  same  definitions  as  before.  Therefore, 
if  we  define 


■■i  = I 'Wk>“ 

k=l 


^2 

r.  = (M  SL  B )“ 
2 m n,m 


(44) 


= r + G“V(i-p)  - x*r  X* 

3X  

Again  SeFy,  and  the  remarks  given  for  Equation  (28)  apply. 


(45) 


The  general  forms  given  above  used  the  symbol  "S"  to  indicate  a 
matrix  function  in  the  resultant  expression.  Furthermore,  the  form  of 
the  result  was  restricted  among  several  possibilities  so  that  SeFm.  This 
restriction  assures  the  uniqueness  of  the  result.  More  formally,  S can 
be  defined  as  follows.  Let  QsFm,  q=trQ,  then 


(46) 
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For  example,  if  Q=e  , then 
dQ  = If  dY 


and 


tr  dQ  = tr[d(e^)]  = tr[d(I+Y+  ~ + 

^ • 


•)] 


= tr[(I+Y  + ^ +'*")dYj  = trCe^dY] 


(47) 


and  3q/3Y=eY.  Clearly,  comparable  results  follow  for  each  element  in  Fm. 

Throughout  the  development  which  resulted  in  the  general  equation  (45), 
attention  has  been  confined  to  the  functions  in  Fm*  It  is  clear  that  any 
function  expressible  as  a convergent  series  has  a matrix  counterpart  [7, 
theorem  4,  p 46].  Thus,  the  set  Fm  could  be  expanded  considerably;  it 
would  include,  among  others,  such  functions  as  logY,  Bessel  functions, 
etc.,  in  which  S can  be  substituted  as  shown  in  Equations  (44),  based 
on  the  result  demonstrated  by  (47),  SeFm» 


Based  on  (47),  the  derivative  of  f(Y)  with  respect  to  Y is  the  same 
as  if  Y were  a scalar.  Thus,  if  Q=sinY, 


gf  = jY  tr(sinY)  = cosY, 


etc. 


Of  particular  importance  in  applications  are  the  matrix  functions 
AXB  and  AX*B.  Table  1,  below  lists  the  gradient  of  the  trace,  3q/3XT  of 
these  two  forms  for  the  functions  defined  in  Fm*  Capital  letters  repre- 
sent matrices  in  Sa,  with  XeSx  and  X*eS*. 
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TABLE  1 


1 


Y 

Q 

3q/3XT 

Y=ABX 

Y 

BA 

CY 

BCA 

Be^A 

sinY 

B(cosY)A 

cosY 

-B(sinY)A 

sinhY 

B(coshY)A 

coshY 

B(sinhY)A 

Y=AX*B^ 

Y 

g"^aV(i-p)-x*bax* 

CY 

g"UVb^(i-p)-x*bcax* 

Y 

e 

G"^  A^ (e^ )^B^ ( I -P ) -X*Be^AX* 

SinY 

G"’a^(cosY^)B  (I-P)-X*B(cosY)AX* 

cosY 

-G"^A^(siny''’)B^Cl-P)»X*B(sinY)AX* 

sinhY 

G'^A^(coshY^)B^(I-P)-X*B(coshY)AX* 

coshY 

G'’A^(sinhY^)B^(I-P)-X*B(sinhY)AX* 

^ G=(XTX),  p=xx* 


More  generality  is  available  using  Equations  (26),  (44),  and  (45). 
Table  2 lists  the  functions  QeF(^,  and  S=3q/3Y  which  are  required  in  the 
gradient  of  the  trace  of  Q,  3q/3x'^,  for  each  choice  of  Q.  Table  3 lists 
a number  of  forms  for  Y and  3q/3XT  in  which  the  required  function,  S,  is 
found  in  Table  2 depending  on  Q. 

TABLE  2 


Q 

AY 

eY 

cosY 

sinY 

coshY 

SinhY 

S 

A 

eY 

-sinY 

cosY 

sinhY 

coshY 
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TABLE  3 


Y 

3q/3XT 

AXB 

BSA 

AXBXC 

BXCSA  + CSAXB 

AXBX^C 

BX^CSA  + (CSAXB 

XX^ 

X^(S+S^) 

XBXJ 

Bl^S  + B^X^S^  j 

AXX  B 

x'bSA  + (BSAX) 

AXXB 

XBSA  + BSAX 

AX*B 

G'^(BSA)^(I-P)  - X*BSAX* 

AX*BX*C 

G"V^(I-P)-X*rX*,  where 
r = BX*CSA  + CSAX*B 

AX*BX*^C 

= G'^r^(I-P)  - X*rX*,  where 
r = BX*^CSA  + (CSAX*B)^ 

= G‘’(S+S^)X*(I-P)  - X*X*^(S+S^)X* 

X*X*^ 

X*X* 

= G"^(sV+X*V)  - X*(X*S+SX*)X* 

X*^X*^ 

= G’^(X*VsX*^)  - X*(sV+X*S^)X* 

= CSAX*B  + G’V^(I-P)  - X*rX*.  where 
r = BXCSA 

AX*BXC 

AX*BX^C 

= (CSAX*B)''’  + G"^r^(I-P)  - X*rX*.  where 
r = BX^CSA 

AX*^BXC 

= CSAX*^B  + G"’r^(I-P)  - x*rX*.  where 
r =•  (BXCSA 
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IV.  APPLICATIONS  TO  KALMAN  FILTERING  THEORY.  This  section  presents  an 
application  of  gradient  techniques  to  the  continuous  Kalman  filter. 

The  optimal  state  vector  defined  by  the  Kalman  filter,  depending  on 
context,  can  be  considered  as  an  estimation  of  the  state  or  as  a filtered 
variable.  Although  there  are  many  ways  one  can  formulate  estimation  or 
filtering  equations,  the  Kalman  formulation  is  sufficiently  general  so 
that  a considerable  variety  of  such  problems  are  expressible  in  that  man- 
ner. The  derivation  of  the  optimal  gain  matrix  is  ideally  suited  to 
demonstrate  the  use  of  the  techniques  developed  in  this  paper. 

The  presentation,  below,  omits  most  details  of  the  theory  already 
extensively  recorded  in  the  open  literature;  sufficient  detail  is  pre-  *' 
sented  so  as  to  obtain  the  matrix  Riccati  equation  for  the  propagation  of 


the  covariance  of  the  state. 

Consider  a "truth  model" 

x(t)(iN  = A x(t)(nS.  + f(t)(i^  t B u(t)(^  (1) 

for  a countably  infinite  sequence  of  continuous  trajectories,  jal, •••,“. 
The  measurement  vector  is 

2(t)(n^.  = A (t)  x(t)(i^.  + v(t)(i^.  (2) 

^ mxn  ^ ^ 

The  structure  of  the  optimal  state  estimate  is  , 

A(t,t)^^  = Ai?(t,t)^^  + f(t^  + W(t)  x(t,t^j  (3) 

where  the  error  in  observation  is 

z(t,t)  ^ = z(t)  ^ - z(t,t)  ^ (4). 

The  error  dynamics  are 

x(t,t^j  = (A-WH)  x(t,t^^  + Bi^^  + (5) 

The  variance  in  the  state  estimation  is 

P(t,t)  » E(x(t,t^^x(t,t))  (6) 
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'1 


with  dynamics  given  by 

^P  ^P 

P = (A-WH)P  + P(A  -H  W ) + BQB  + WRW  (7) 

where  the  standard  uncorrelated  assumptions  are  assjmed  and 
E = Q(t,x)  6(t,t) 

and 

E ” R(t,x)  6(t,x) 

Any  full  rank  W used  in  Equation  (3)  will  generate  a trajectory;  the 
problem  is  to  find  the  W(t)  that  will  minimize  the  trace  of  Equation  (7). 
It  can  be  shown  that  if  P is  minimized,  then  P is  minimized,  hence 


3P  . 3 


(-WHP) 


3W  3W 

From  Table  3,  Section  III 


-Ajr  (PhV)  + -Ajr  (WRW'^) 

3W  3W 


(8) 


3W 

3 


3W 

3 


3W 


tr  (-WHP)  = - HP 


tr  (PhV)  = + HP 


tr  (WRw'^)  = 2Rw'^ 


(9) 

(10) 
(11) 


Using  Equations  (9)  through  (11)  in  Equation  (9),  gives 


3P  T 

- 0 = - 2HP  + 2RW 


3W 

from  which 

w’^  = r"^p 
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That  is , 


T -1 

W = PH  R 


(12) 


Which  is  the  standard  Kalman  gain  matrix. 
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APPENDIX 


NOTATION 


1.  Upper  case  Roman  letters  A,  B,  C,  •••,  generally  refer  to  rectan- 
gular matrices. 


Special  Forms: 


I 


Is  a unit  matrix; 


I 

P 


a pxp  unit  matrix. 


0 


Is  a matrix  of  zeros; 


0 

P 


a pxp  matrix  of 


zeros . 


X A rectangular  matrix  of  random  variables  or  a set  of  mea- 

surements. Each  column  is  a measurement  and  the  number  of 
rows  is  generally  greater  than  or  equal  to  the  number  of 
columns.  X is  norroally  "full-rank" , that  is,  the  rank  of  X 
equals  the  number  of  columns . 


X*  Is  the  generalized  inverse  of  X.  In  particular,  it  is  the 

pseudo  inverse  of  X when  X is  assumed  to  be  full-rank . 

That  is,  if  the  rank  of  X is  equal  to  the  number  of  its 
columns,  then  X*=(xTx)"lxT.  (NOTE:  If  the  rank  of  X were 

equal  to  the  numbers  of  its  rows,  then  X*»X^(XXT)-1 . ) If 
X is  square,  then  X*=X-1. 

. . T 

G Is  the  Grammian  of  a full  row  rank  matrix;  thus  G=X  X* 


W Generally  refers  to  a matrix  of  weights. 


P,Q,R  Frequently  refer  to  the  covariances  of  the  optimal  estimate 

of  the  state-vector,  process  noise,  and  measurement  noise,  resp. 

P,P  Represent  a projector  of  a variable  array  into  its  range 

space  and  null  space,  respectively.  Both  P and  P are 

idempotent  and  Pp=0. 

2.  Lower  case  Roman  letters,  a,  b,  c,  •••,  refer  to  vectors. 

Special  Forms.  Frequently  an  upper  case  Roman  letter,  say  A,  will 
refer  to  an  unspecified  array.  For  those  special  cases  in  which  the 

array  is  to  be  a vector,  lower  case  a will  replace  A. 
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3.  Operators. 

Superscripts 

T 

T Is  the  transpose  of  an  array;  thus  A . 

-1  Is  the  inverse  of  a square  array;  thus  A 

Is  the  pseudo  (generalized)  inverse  of  an  array;  thus  A*. 
(If  A has  full -rank,  A*=(A^A)  ^A^.) 


Symbols 

tr 

< •> 


bq, 


Is  the  trace  of  a square  array;  thus  trA=  ^ a.,  where 
A=(a^^)  is  an  matrix.  i=l 

The  Dirac  bra  and  ket,  respectively,  represent  row  and 
column  vectors  respectively.  Since,  in  this  paper,  all 
arrays  are  of  real  numbers , only , the  inner  product  of  two 
vectors  a and  b is  <ab>=<ba>  or  <aMb>=<bMTa>  over  the  metric 
array , M . 

Represent  row  and  column  vectors  a and  b,  having  p and  q 
compontents,  respectively. 


^ a,6  natural  numbers,  are  used  to  represent  arrays  A and  B 

partitioned  into  their  column  space  or  row  space, respectively, 
or  are  used  to  emphasize  that  two  vectors,  a and  b,  are 
defined  in  a "space"  and  its  "cospace"  defined  by  a "basis" 
and  "reciprocal  basis". 

Appearing  above  an  array  respresent  an  optimal  estimate  and 
its  optimal  error,  respectively,  thus  X and  X.  Also  note 
that  such  an  array  can  be  written  as  X-X+X. 

E(-),  Represent  expected  value  and  covariance  of  the  arguments, 
cov(‘)  respectively. 
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